// 版权归Go作者所有。版权所有。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

/*
	Bessel function of the first and second kinds of order n.
*/

// 原始C代码和下面的长注释是
// 来自FreeBSD的/usr/src/lib/msun/src/e_jn.C和
// 随此通知而来。go代码是原始C.
// 
// ===================================================================================================
// /Sun Microsystems，Inc.1993年版权所有。
// 
// 由SunPro开发，SunPro是Sun Microsystems，Inc.的一家企业。
// 使用、复制、修改和分发本
// 软件的权限是免费授予的，前提是保留本通知
// 的权利。=========================================================================================================
// /
// ieee754 jn（n，x），ieee754 yn（n，x）
// /第一类和第二类浮点贝塞尔函数
// /第n阶阶
// y0（-ve）=y1（-ve）=yn（n，-ve）是带有无效信号的NaN。
// 注2.关于jn（n，x），yn（n，x）
// 对于n=0，调用j0（x），对于n=1，j1（x）对于n<x，使用前向递归从j0（x）和j1（x）的值开始，使用前向递归。对于n>x，使用连分式近似值，表示
// j（n，x）/j（n-1，x）对j（n，x）进行求值，然后从假设值
// 递归。j（0，x）的结果值为
// 开始使用反向
// 与实际值进行比较，以更正
// j（n，x）的假设值。
// 
// 所有
// n>1的值都使用前向递归外，在所有方面都是类似的。

// Jn返回第一类的n阶贝塞尔函数。
// 
// 特殊情况是：
// Jn（n，±Inf）=0 
// Jn（n，NaN）=NaN 
func Jn(n int, x float64) float64 {
	const (
		TwoM29 = 1.0 / (1 << 29) // 2**-29 0x3e10000000000000 
		Two302 = 1 << 302        // 2**302 0x52D0000000000000 
	)
	// 特殊情况
	switch {
	case IsNaN(x):
		return x
	case IsInf(x, 0):
		return 0
	}
	// J（-n，x）=（-1）**n*J（n，x），J（n，-x）=n*J（n，x）
	// 因此，J（-n，x）=J（n，-x）

	if n == 0 {
		return J0(x)
	}
	if x == 0 {
		return 0
	}
	if n < 0 {
		n, x = -n, -x
	}
	if n == 1 {
		return J1(x)
	}
	sign := false
	if x < 0 {
		x = -x
		if n&1 == 1 {
			sign = true // 奇数n和负数x 
		}
	}
	var b float64
	if float64(n) <= x {
		if x >= Two302 { // x>2**302 

			// （x>>n**2）
			// /Jn（x）=cos（x-（2n+1）*pi/4）*sqrt（2/x*pi）
			// 
			// n sin（xn）*sqt2 cos（xn）*sqt2 
			// /------------------------------------
			// 0 s-c c+s 
			// /1-s-c-c+s 
			// 2-s+c-c-s 
			// 3 s+c-s 

			var temp float64
			switch s, c := Sincos(x); n & 3 {
			case 0:
				temp = c + s
			case 1:
				temp = -c + s
			case 2:
				temp = -c - s
			case 3:
				temp = c - s
			}
			b = (1 / SqrtPi) * temp / Sqrt(x)
		} else {
			b = J1(x)
			for i, a := 1, J0(x); i < n; i++ {
				a, b = b, b*(float64(i+i)/x)-a // 避免下溢
			}
		}
	} else {
		if x < TwoM29 { // x<2**-29 
			// x很小，返回J（n，x）
			// J（n，x）=1/n！*（x/2）**n-…

			if n > 33 { // 下溢
				b = 0
			} else {
				temp := x * 0.5
				b = temp
				a := 1.0
				for i := 2; i <= n; i++ {
					a *= float64(i) // a=n！
				}
				b /= a
			}
		} else {
			// 使用后向递归
			// x**2x**2abcdefg
			// 2n-2（n+1）-2（n+2）
			// 
			// 1 
			// （对于大x）=---------------…
			// 2n2（n+1）2（n+2）如果w=2n/x，h=2/x，则上述商等于连分数：
			// 1 
			// 1 
			// w----
			// 1 
			// w+h----
			// w+2h-…
			// 
			// 确定需要多少个术语，让
			// Q（0）=w，Q（1）=w（w+h）-1，
			// Q（k）=（w+k*h）*Q（k-1）-Q（k-2），
			// Q（k）>1e4适用于单
			// Q（k）>1e9适用于双
			// Q（k）>1e17适用于四

			// /确定k 
			w := float64(n+n) / x
			h := 2 / x
			q0 := w
			z := w + h
			q1 := w*z - 1
			k := 1
			for q1 < 1e9 {
				k++
				z += h
				q0, q1 = q1, z*q1-q0
			}
			m := n + n
			t := 0.0
			for i := 2 * (n + k); i >= m; i -= 2 {
				t = 1 / (float64(i)/x - t)
			}
			a := t
			b = 1
			// 因此，如果n*（log（2n/x））>。。。
			// 单体8.8722839355e+01 
			// 双体7.09782712893383973096e+02 
			// 长双体1.13565234062949491931077970765006170E+04 
			// 可能下溢到零

			tmp := float64(n)
			v := 2 / x
			tmp = tmp * Log(Abs(v*tmp))
			if tmp < 7.09782712893383973096e+02 {
				for i := n - 1; i > 0; i-- {
					di := float64(i + i)
					a, b = b, b*di/x-a
				}
			} else {
				for i := n - 1; i > 0; i-- {
					di := float64(i + i)
					a, b = b, b*di/x-a
					// 避免标度b伪溢出
					if b > 1e100 {
						a /= b
						t /= b
						b = 1
					}
				}
			}
			b = t * J0(x) / b
		}
	}
	if sign {
		return -b
	}
	return b
}

// Yn返回第二类的n阶贝塞尔函数。
// 
// 特殊情况为：
// Yn（n，+Inf）=0 
// Yn（n≥ 0，0）=-Inf-
// Yn（n<0，0）=+Inf如果n是奇数，-Inf如果n是偶数
// Yn（n，x<0）=NaN 
// Yn（n，NaN）=NaN 
func Yn(n int, x float64) float64 {
	const Two302 = 1 << 302 // 2**302 0x52D0000000000000 
	// 特殊情况
	switch {
	case x < 0 || IsNaN(x):
		return NaN()
	case IsInf(x, 1):
		return 0
	}

	if n == 0 {
		return Y0(x)
	}
	if x == 0 {
		if n < 0 && n&1 == 1 {
			return Inf(1)
		}
		return Inf(-1)
	}
	sign := false
	if n < 0 {
		n = -n
		if n&1 == 1 {
			sign = true // 如果n<0&&&n | n |奇数
		}
	}
	if n == 1 {
		if sign {
			return -Y1(x)
		}
		return Y1(x)
	}
	var b float64
	if x >= Two302 { // x>2**302 
		// ------------------------------------
		// 0 s-c-c+s 
		// /1-s-c-c+s 
		// /2-s+c-c-s 

		var temp float64
		switch s, c := Sincos(x); n & 3 {
		case 0:
			temp = s - c
		case 1:
			temp = -s - c
		case 2:
			temp = -s + c
		case 3:
			temp = s + c
		}
		b = (1 / SqrtPi) * temp / Sqrt(x)
	} else {
		a := Y0(x)
		b = Y1(x)
		// /如果b是-inf 
		for i := 1; i < n && !IsInf(b, -1); i++ {
			a, b = b, (float64(i+i)/x)*b-a
		}
	}
	if sign {
		return -b
	}
	return b
}
